Normalisation
NadineBestard
08/03/2021
Set-up
library(here) # for reproducible paths
library(SingleCellExperiment)
library(scater) # For qc
library(scran) # For normalisation
library(Matrix) # For log transorming the raw data
library(ggplot2) # To add titles to plots
library(batchelor) # Mnn, for batch correction# Adapted function from VISION to log tranform sparse matrix
# I could not download the package
matLog2 <- function(spmat, scale = FALSE, scaleFactor = 1e6) {
if (scale == TRUE) {
spmat <- t( t(spmat) / colSums(spmat)) * scaleFactor
}
if (is(spmat, "sparseMatrix")) {
matsum <- summary(spmat)
logx <- log2(matsum$x + 1)
logmat <- sparseMatrix(i = matsum$i, j = matsum$j,
x = logx, dims = dim(spmat),
dimnames = dimnames(spmat))
} else {
logmat <- log2(spmat + 1)
}
return(logmat)
}Normalisation by deconvolution
In order to correct for systematic differences in sequencing coverage between libraries we will normalise the dataset. This involves dividing all counts for each cell by a cell-specific scaling factor, often called a “size factor” (Anders and Huber 2010). The assumption here is that any cell-specific bias (e.g., in capture or amplification efficiency) affects all genes equally via scaling of the expected mean count for that cell. The size factor for each cell represents the estimate of the relative bias in that cell, so division of its counts by its size factor should remove that bias.
Specifically we will used the deconvolution method available in the scran package. This method allows to take in consideration the composition bias between samples (Lun et al., 2016)
# Only compute if first time
if (!(file.exists(here("processed", "sce_norm.RDS")))) {
sce <- readRDS(here("processed", "sce_QC.RDS"))
# For reproducibility
set.seed(100)
# Quick clustering to pool samples together and deal with 0 counts
quick_clusters <- quickCluster(sce)
# Calculate size factors
sce <- computeSumFactors(sce, cluster = quick_clusters, min.mean = 0.1)
# Check that there are not negative size factors
summary(sizeFactors(sce))
# Apply size factors and log transform them
sce <- logNormCounts(sce)
# Also log normalise the raw counts
assay(sce, "logcounts_raw") <- matLog2(counts(sce))
saveRDS(sce, here("processed", "sce_norm.RDS"))
} else{
sce <- readRDS(here("processed", "sce_norm.RDS"))
}On top of normalisation the data it is also log-transformed. The log-transformation is useful as differences in the log-values represent log-fold changes in expression. Or in other words, which is more interesting - a gene that is expressed at an average count of 50 in cell type A and 10 in cell type B, or a gene that is expressed at an average count of 1100 in A and 1000 in B? Log-transformation focuses on the former by promoting contributions from genes with strong relative differences.
Assess Confunding factors impact
Variance Explained plots
Variable-level metrics are computed by the getVarianceExplained() function (before and after normalization). This calculates the percentage of variance of each gene’s expression that is explained by each variable in the colData of the SingleCellExperiment object. We can then use this to determine which experimental factors are contributing most to the variance in expression. This is useful for diagnosing batch effects or to quickly verify that a treatment has an effect.
The percentage of variance explained by a factor is on the x axis, and in the y axis there is the density of the R-squared values across all genes. The “total” is the total number of molecules, that correlates with the detected number of genes, detected.
Before normalisation
Before normalisation it is expected that most variance will be explained by the sequencing depth, i.e. the total numbber of umis and the total number of genes
# Before normalisation
# Only compute if first time
if (!(file.exists(here("processed", "variance_explained.RDS")))) {
# Calculate the matrix (time consuming step)
var <- getVarianceExplained(
sce,
exprs_values = "logcounts_raw",
variables = c(
"tissue",
"chip",
"Sample",
"age",
"subsets_mt_percent",
"detected",
"total"
)
)
saveRDS(var, here("processed", "variance_explained.RDS"))
#If not just load created object
} else{
var <- readRDS(here("processed", "variance_explained.RDS"))
}
plotExplanatoryVariables(var)## Warning: Removed 126 rows containing non-finite values (stat_density).
After normalisation
We can see how there is less variance explained now by factors such as the detected genes or the number of counts
# After normalisation
if (!(file.exists(
here("processed", "variance_explained_norm.RDS")
))) {
var_norm <- getVarianceExplained(
sce,
variables = c(
"tissue",
"chip",
"Sample",
"age",
"subsets_mt_percent",
"detected",
"total"
)
)
saveRDS(var_norm, here("processed", "variance_explained_norm.RDS"))
} else{
var_norm <- readRDS(here("processed", "variance_explained_norm.RDS"))
}
plotExplanatoryVariables(var_norm)## Warning: Removed 126 rows containing non-finite values (stat_density).
# Again after normalisation but other parametres
if (!(file.exists(
here("processed", "variance_explained_norm_2.RDS")
))) {
var_norm_2 <- getVarianceExplained(
sce,
variables = c(
"tissue",
"chip",
"age",
"genotype",
"mouse",
"total"
)
)
saveRDS(var_norm_2, here("processed", "variance_explained_norm_2.RDS"))
} else{
var_norm_2 <- readRDS(here("processed", "variance_explained_norm_2.RDS"))
}
plotExplanatoryVariables(var_norm_2)## Warning: Removed 108 rows containing non-finite values (stat_density).
It is reassuring to see that more variance is explained by the mice the sample came from than from the chip.
Dimensional reduction
Another way to assess the variance is by a plot with a PCA plot. Here again we can see how the sequencing depth(sum) explains most of the variance before the normalisation (49%)
raw <- runPCA(sce, exprs_values = "logcounts_raw")
plotPCA(raw, colour_by= "chip", size_by="sum") + ggtitle("Before normalisation")
sce <- runPCA(sce)
plotPCA(sce, colour_by= "chip", size_by="sum") + ggtitle("After normalisation")
plotPCA(sce, colour_by= "chip", point_size=0.1) +
ggtitle("After normalisation, small dots")
Batch correction with MNN
Combat was shown to outperform other batch correction methods for simple batch correction (Buttner et. al, 2019). However, this will also regress other biological differences that are not well balanced between batches. Integration techniques account for this fact, with the downside that it can lead to over-correction due to increased degrees of freedom of these non-linear methods.
We use a merge tree, useful for merging together batches that are known to be more closely related before attempting difficult merges involving more dissimilar batches.
set.seed(100)
sce_batch_corrected <- fastMNN(sce,
batch = factor(sce$chip),
merge.order = list(
list(list("3","5"), list("4","6")),
list("1","2")
)
)
# Todo, try regressBatches() with the design fromcolData(sce_batch_corrected) <- cbind(colData(sce), colData(sce_batch_corrected))
sce_batch_corrected <- runPCA(sce_batch_corrected, exprs_values = "reconstructed" )
plotPCA(sce_batch_corrected, colour_by= "chip", point_size=0.1) +
ggtitle("After batch correction small dots")
Session Info
Click to expand
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
## [5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
## [7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] batchelor_1.6.2 Matrix_1.3-2
## [3] scran_1.18.5 scater_1.18.6
## [5] ggplot2_3.3.3 SingleCellExperiment_1.12.0
## [7] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [9] GenomicRanges_1.42.0 GenomeInfoDb_1.26.2
## [11] IRanges_2.24.1 S4Vectors_0.28.1
## [13] BiocGenerics_0.36.0 MatrixGenerics_1.2.1
## [15] matrixStats_0.58.0 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] locfit_1.5-9.4 Rcpp_1.0.6
## [3] rsvd_1.0.3 lattice_0.20-41
## [5] rprojroot_2.0.2 digest_0.6.27
## [7] utf8_1.1.4 R6_2.5.0
## [9] evaluate_0.14 highr_0.8
## [11] bluster_1.0.0 pillar_1.5.1
## [13] sparseMatrixStats_1.2.1 zlibbioc_1.36.0
## [15] rlang_0.4.10 irlba_2.3.3
## [17] rmarkdown_2.7 labeling_0.4.2
## [19] BiocNeighbors_1.8.2 statmod_1.4.35
## [21] BiocParallel_1.24.1 stringr_1.4.0
## [23] igraph_1.2.6 RCurl_1.98-1.2
## [25] munsell_0.5.0 beachmat_2.6.4
## [27] DelayedArray_0.16.2 compiler_4.0.4
## [29] vipor_0.4.5 BiocSingular_1.6.0
## [31] xfun_0.21 pkgconfig_2.0.3
## [33] ggbeeswarm_0.6.0 htmltools_0.5.1.1
## [35] tibble_3.1.0 gridExtra_2.3
## [37] GenomeInfoDbData_1.2.4 edgeR_3.32.1
## [39] fansi_0.4.2 viridisLite_0.3.0
## [41] crayon_1.4.1 withr_2.4.1
## [43] bitops_1.0-6 grid_4.0.4
## [45] gtable_0.3.0 lifecycle_1.0.0
## [47] magrittr_2.0.1 dqrng_0.2.1
## [49] scales_1.1.1 stringi_1.5.3
## [51] scuttle_1.0.4 farver_2.1.0
## [53] XVector_0.30.0 viridis_0.5.1
## [55] limma_3.46.0 DelayedMatrixStats_1.12.3
## [57] ellipsis_0.3.1 vctrs_0.3.6
## [59] tools_4.0.4 glue_1.4.2
## [61] beeswarm_0.3.1 ResidualMatrix_1.0.0
## [63] yaml_2.2.1 colorspace_2.0-0
## [65] knitr_1.31